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ABSTRACT 


This report summarizes the analyses of deflagration to detonation 
transition [DDT) occurring in a packed bed of granular, high energy solid 
propellant. A reactive two-phase flow model of this phenomena is solved 
by utilizing a lax-Wendroff finite differencing technique. A brief 
overview of the well known shock jump conditions for one-dimensional, 
one-phase flow with heat addition is reported, and a similar analysis 
for one-dimensional, two-phase reactive flow is discussed. Improvements 
made in the gas phase nonideal equation of state, gas permeability, and 
numerical integration t chniques allow for the prediction of a transition 
to a steady detonation from initiation by deflagration. 

Analyses are presented that clearly indicate the effect of the 
propellant physical and chemical parameters on the predicted run-up 
length to detonation. Predictions of this run-up length to detonation 
are presented as a function of propellant chemical energy, burning rate, 
bed porosity, and granulation (size). Limited comparison with actual DD1 
data in the literature indicates good qualitative agreement with these 
predictions. 
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CHAPTER ONE 


DEF LAGRATION-^0-DETONATION TRANSITION 

1.1 Introduction 

This report summarizes the analysis associated with the accelerating 
deflagration wave in a porous medium of reactive solid propellant. The 
phenomenon of DDT (.deflagration-to-detonation transition) in solid pro¬ 
pellants, e pecially solid propellants burning in rocket motor environ¬ 
ments, is not usually considered a hazard. However, it. may be that un¬ 
der certain situations, for example a grain structure failure, the solid 
motor may crack and form regions of granular or porous propellant. When 
flame from the surface deflagrating propellant reaches this seam of porous 
material it will accelerate into this medium and be supported by convective 
heat transfer from the burnt gas into the unignited porous region. If in 
addition to this the product gases are confined to a finite volume, the 
accelerating deflagration could transit into a detonation. Propellants 
exhibiting a high chemical energy per unit mass and capable of rapidly 
generating gases through their burning rate are more likely to experience 
this type of deflagration to detonation transition. 

Analysis of the flow.* in such an unsteady two-phase mixture is a com¬ 
plicated exercise. Work has been underway at the University of Illinois 
since 197S to develop a reactive hydrodynamic code in which the combus¬ 
tion of porous propellants can be modeled in such a way as to predict the 
behavior of a convectively driven flame in a confined situation. Details 
of these modeling exercises are found in References 1-4. In an evolutionary 
mariner this work included the formulation of the two-phase flow conserva¬ 
tion equations, first assuming that the mixture was a continuum, and at a 
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later date treating each phase a:, a separate continuum irrespective oz 
the mixture properties. The most recent analysis of the unsteady two-phase 
flows associated with DDT is documented in the paper by Hoffman and Krier 
[3]. This reference, therefore, represents the smarting point for the 
work that is presented here. 

The reader, after reviewing the above noted reference, will understand 
that the modeling of this transient phenomena utilizes a number of impor¬ 
tant constitutive relations, which form closure of the conservation equa¬ 
tions. For example, one must have information on the burning rate of the 
material that is a function of the surrounding pressure and temperature. 

One must also have laws for the dynamic gas permeability and the subse¬ 
quent heat transfer rates of the hot gases as they are forced into the 
unignited porous material. In addition, relations which represent the 
resistance to compaction of the solid matrix must be included. Equations 
of state, not only for the high pressure in the product gases, but also 
for the solid itself, must be supplied (as shown in the work by Hoffman 
and Krier). The assumption of an incompressible solid, although provid¬ 
ing some reasonable answers s far as the deflagration speed, is not an 
accurate indicator of the peak pressures that are possible during the 
accelerating deflagration mode. Since these pressures are precursor to 
the final detonation solutions that would be expected, it is clear that 
a compressible solid must be modeled. 

In summary, the analysis of the DDT problem requires the solution of 
the conservation equations in both phases and the necessary constitutive 
relations. The conservation equations form a system of nonlinear hyper¬ 
bolic equations, which require numerical finite differencing schemes. 

This report will summarize several of these integration schemes and will 
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evaluate which are more useful for this reactive flow problem. There are, 
of course, many numerical techniques available to handle hyperbolic par¬ 
tial differential equations which allow for solutions of shock waves. 
Obviously, not all of these have been treated in this work. 

1.2 Previous Results on DDT Modeling 

A review of References 2-4 indicates that a steady state detonation 
solution was not a predicted result. As will be shown in Sections 1.3 
and 2.4, for propellant chemical energy, the initial solids loading con¬ 
sidered and the material parameters (e.g. y) , a steady state detonation 
(CJ) would propagate at speeds of the order of 5-8 nun/ysec, with a deto¬ 
nation pressure of the order of 12-20 GPa. Although a fairly rapid flame 
front, often approaching a steady state speed of 2 mm/ysec was typical of 
the solutions presented in the work by Gokhale and Krier [2], Kezerle and 
Krier [4], and Hoffman and Krier [3], the associated peak pressures were 
never of significant manner to suggest that these final flame spreading 
rates were characteristic of an actual detonation, since fronts ranging 
from j to 2 mm/ysec cannot support a detonation (see following section). 
However, the peak pressures previously reported, which ranged between 
1/2 to 5 GPa, were consistent with the velocities. 

Continued work, which is reported here, indicates that it is now 
possible to obtain a detonation solution, but in order to do so a num¬ 
ber of important modifications and corrections are required. To under¬ 
stand what these changes are and why they are necessary, it is appro¬ 
priate at this point to first provide a review on the topic of detona¬ 
tion. Although the discussion that follows is documented in various 
manners in several textbooks, including such a review here will allow 
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for a clearer understanding of the steady state solution that is termed 
"detonation", as well as define the flow properties of the DDT analysis. 

1.3 Jump Conditions for Reactive Flow 

In order to understand the jump conditions across a detonation wave 
in two-phase (sol?d-gas) reactive flow, one should first look at the solu¬ 
tion of the one-phase (all gas] flow problem. Pickett and Davis [5], anu 
Strehlow [6] -'evelop these jump equations in more detail in their texts. 
This section will give the reader a brief review of their work and high¬ 
light some of the key assumptions made in developing the Hugoniot and 
Rayleigh line equations. 

In Figure 1 a detonation wave with velocity D is moving through a 
constant area duct into the unignited or "cold" end which is at rest. 

For this analysis all chem cal reactions are assumed to take place in a 
narrow reaction zone. The products of combustion are shown moving with 
velocity u in the same direction as D. In this figure, the subscript A 
on pressure, density, and velocity designates the unignited or reactants 
side, and B indicates the product side of the combustion zone. 

To better understand the jump conditions across the shock in Figure 
1, a new frame of reference may be helpful. From a coordinate system 
attached to the moving detonation wave, an observer located at the origin 
would see the reactants moving with velocity u^ = D into the stationary 
reaction zone and the products of combustion exiting with velocity 
Ug => (D-u) (see Figure 2). 

Conservation of ass and momentum through the flow area shown in 
Figure 2 gives respectively 

P A U " P B U B = P B (D “ U) 


( 1 ) 








5 


Moving Reaction Zone 7 



Fig. 1. Detonation Wave Moving into Stationary Reactants. 
(Reference Frame Fixed). 
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and 


Pa ° 2 ♦ 


P A 


P B U B 


+ P B - 


P b (D-uJ 


+ P, 


( 2 ) 


In this model all viscous effects are confined within the bounds of the 
shock discontinuity. The expression for conservation of momentum (Eq. 2) 
can be simplified by making use of Equation 1 to obtain 


P 


B 


P A ' V” 


( 5 ) 


Later, it. will be assumed that P A is nagligible with respect to P^, but 
for the time P A will be retained in the momentum equation. 

Elimination of u, the reaction products velocity expressed in the 
fixed coordinate system, from Eqs. 1 and 3 gives the well known Rayleigh 
or Michelson line. 

p/d" - (P D - P a )/(v a - v B ) * 0 (4) 


Figure 3 illustrates a Rayleigh line on a P-v diagram for one specific 

2 2 

value of D. The line has slope -p A D and passes through the point C? A » V A )• 

For the case where D -*■ the Rayleigh line approaches vertical and for 

D -*■ 0, the Rayleigh line is horizontal. In order to satisfy the conserva¬ 
tion of mass and momentum, the final state fP . v 1 must also lie nn the 

.' ' ' “ B' ' B' •” 

Rayleigh line as shown in Figure 3. 

Simultaneously solving Equations 1 and 3 by eliminating D instead of 
u, the result would be an expression for u, the gaseous product velocity 
in the fixed reference frame. 


uZ * (P B - P A ) (V A - V (*) 

Figure 4 shows several curves of constant u plotted on a P-v diagram. 
Solving Equations 1 and 2 for = (D-u) one would obtain 



7 



Fig, 3. Illustration of Rayleigh lir.e. 
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V B 2 ( P B - P ^‘ V A 


V 


(5a) 


Equation 5a will be important later in showing that for a CJ detonation 
the product velocity is sonic with respect to the detonation front. 

The third conservation equation, the conservation of energy across 
the reaction zone, can be written 


E. „ + P.v. + 4 D 2 * E. „ + P.v D + i u D 2 = E. „ 

int. A A 2 int Q B B 2 B mt n 

A B B 


* P B V B * I (D ' U)2 


( 6 ) 


In Equation 6, represents the specific internal energy at each re¬ 
spective location and v is the specific volume v = 1/p. Since C v> the 
specific heat at a constant volume, is assumed constant throughout the 
process, the internal energy at location B is simply 


E. = C T_ 
v 8 


(7a) 


Likewise, the internal energy at location A is similar with the inclusion 
of the chemical energy yet to be released. 


E. = C T + E 
int. v A CHF.M 
A 


(7b) 


By convention a positive value of ^ is endothermic and a negative value 


ovorhormi< 


Elimination of u and D from the energy equation (Eq. 6) is obtained by 
substituting in the mass and momentum equations (I'qs. 1 and 3). The re¬ 
sult is a relation between t. P, and v in both reactant and prcduct 

int’ r 

states known as the Hugoniot equation. It is expressed as 



t int 2 (P B 
A 


V ( V A - 


V = 


( 8 ) 


If one assumes that is a constant and the reaction products obey an ideal 
equation of state 
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Pv = RT 

then the internal energy expressions (Eqs. 7a and 7b) take the form 


(9) 


E. , 
mt r 


C P v P v 

v B B B B 

R = (Y-l) 


and 


C P,v. 
r- V A A. 

E. = —-- + E 

int. R 
A 


= VA 

CHEM (y-l) + h CHEM 


(10a) 


(10b) 


where y is the ideal gas limit of Substitution of Equations 10a 

and 10b into Equation 8 gives a new form of the Hugoniot equation. 


P 

+ izi 


- 1 + 

Izl 

2 

a 

Xii 

rH 

+ 

< 

cu 

-1 

l V A H 


Y + 1 


Y+l 

L 1 


2E 


CHEM 


P. V A 
A A 


( 11 ) 


On a P-v diagram Equation 11 plots as a hyperbola passing through (P , v^) 
for the case where E^^ = 0 and displaced from this point for increasing 
values of E^^ (see Fig. 5). From the solution of the conservation equa¬ 
tions, the Hugoniot relation states that flow coming in at pressure P = P 

A 

and specific volume v = V A with chemical energy E^^ = E CHEM , must have 

a final state (P , v ) which lies on the Hugoniot. As discussed 

o D LrihM^ 

earlier, the Rayleigh line connecting the initial state (P , v ) with the 

A JK 

final state (Pg, v^) dictates a unique solution for the detonation velocity, 
D, for that process. Therefore, the steady s ate flow process between 
states A and B shown in Figure 2 can be uniquely modeled by the Rayleigh 
line and Hugoniot curve on Figure 6, knowing the initial state (P,, v ), 
the final stat< (Pg, Vg) and the chemical energy of the re .ctant. 

There are certain areas on the P-v diagram where a final solution to the 
flow process is .impossible. For instance, the pressure and specific 




Hugomot 
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volume cannot both increase over the shock discontinuity. 

Figure 6 shows a Kugoniot curve and two Rayleigh lines, and P 7 , 
on the same P-v diagram. For Rayleigh line labeled R2 there are two solu¬ 
tions labeled S for strong and W for weak. At the strong solution the 
flow downstream of the shock (u D in Fig. 2) is subsonic with respect to 

D 

the stationary shock, and a pressure disturbance initiated, traveling at 
the local speed of sound in the fluid, will propagate back and overtake 
the shock front. Looking at the same strong solution but in the fixed re¬ 
ference frame (moving detonation wave), a pressure fluctuation behind the 
detonation front (i.e., moving the rear wall in a closed system) will over¬ 
take the front and the front will adjust itself to that change. For the 
weak solution the flow behind the front (u g in Fig. 2) moves away from 
die stationary front at a velocity greater than the local sound velocity. 
Therefore, a pressure disturbance in this case will not be felt by the 
front. Both the strong and weak detonations will be discussed in more 
detail in the next section. 

In the case of the Rayleigh line labeled R^ and the Hugoniot curve, 
there is only one unique solution to the conservation equations. At this 
point the two lines are tangent and the final st. te is identified as the 
CJ (Chapman-Jouguet) point. For this case, the flow behind the front at 
the CJ condition is moving at the local speed of sound with respect to 
the front (i.e., u fi = a ) . 

In order to prove that the flow is sonic with respect to the deto¬ 
nation front,one should first look at the thermodynamic definition of 
sound speed. 






( 12 ) 


















In Equation 12 the subscript s indicates the derivative is evaluated 
at constant entropy. Therefore, for a detonation the local sound velocity 
a_ the CJ point is simply obtained by evaluating the slope of an isen- 
trepe through that point on a P-v plot. 

Since the slope of the Hugoniot and Rayleigh lines are equal at the 
CJ point, it is useful to solve both derivatives and equate them. For 
the Rayleigh line 

<f> ■ - °A d2 ■ - ( P B - V-'Cv A - V < 13 > 

K 

and for the Hugoniot curve 

<f> H • 2 ‘f> n ' tv A - V ‘ < P B - P J ' ‘"A - V t») 

Equating the right hand sides of Eqs. 13 and 14 you obtain 

" - P U5) 

H CJ 

From the thermodynamic relation 


dE =■ Tds - Pdv 


(16) 


you obtain 



(17) 


by evaluating Equation 16 on an isentrope. Equations 15 and 17 imply 
that at the CJ point an isentrope is tangent to both Rayleigh and Hugoniot 
lines. From this 


(~) 

'•ev' 


CJ 


fP B - V 
' V A * V 


(18) 


Therefore, at the CJ point the local sound speed can be evaluated by using 
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liquations 12 and 18 to obtain 


2 

* CJ 



rp - p i 
_ _k’ 

(V A - V 


(19) 


This is exactly the same as Equation 5a, the square of the product ve¬ 
locity for a stationary detonation. Therefore, the product velocity for 
a CJ detonation is sonic with respect to the detonation front. 

From Equation 18, one can define gamma to be <the negative logarith¬ 
mic slope of the isentrope through the CJ point. 


y =-(3ln P/31n v) = 

The analysis of detonation processes presented in this section as¬ 
sumes y is not only a constant value but is the same in both product and 
reactant sides of the detonation wave. For most gases Y is about 1.2 [7], 
In the case of a sc lid (explosive) detonating, Fickett and Davis [5] use 
a value of y = 3.0 and assumes the ideal state equation holds for the 
product gases. From this assumption and estimate of y, one could directly 
apply the analysis formulated in this section to the solid detonation. 

However, for t wo-phase reactive flow the mixture momentum and mix¬ 
ture energy differ from Eqs. 3 and 4 by interactive stress and stress 
work terms. Mso, the gaseous products cannoc be assumed ideal when ex¬ 
perimental work shows detonation pressures of the order 15-20 GPa (2.17-2.9 
Mpsi) [8]. Because of the high detonation pressure a nonideal equation 
of state v,;s implemented. This will be discussed in more detail in Chapter 
2. Also, the jump conditions for a mixture with interactive terms is 
outlined in Appendix A. It is similar to the maverial presented in this 
section with the exception of the stress and stress work terms. 
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1.4 Rear Boundary Condition for Steady State Detonations 


Reference 5 shows that the velocity, D, at which a detonation 
propagates into the reactants is dependent on the velocity at which the 
rear wall, or boundary, travels with respect to the CJ product velocity 
u . The product velocity was discussed in the previous section on the 
Rayleigh line and Figure 4 shows several curves of constant product ve¬ 
locity. The CJ product velocity would be the curve intersecting the CJ 
point in Figure 6. Three different possibilities for the movement of 
the rear wall with respect to the CJ product velocity exist and will be 
discussed in this section. 

The first case to consider is where the rear wall velocity is greater 
than the product velocity at the CJ condition. On the P-v diagram (see 
Fig. 6) this would correspond to the solution labeled S for "strong". 
Here, the flow field following the detonation is constant with the pro¬ 
duct velocity u = u w> Examining Figure 6, one can see that the final 
state pressure, P_, is also greater than that of the corresponding CJ con- 

D 

dition. Because the flow behind the detonation front is subsonic with 
respect to the front, any perturbation in the rear- wall will propagate 
forward to the pressure front. Figure 7a shows a pressure profile for an 
overdriven detonation with the CJ pressure labeled on the figure. 

The eff' cts on the pressure by reducing the rear wall velocity can 
be seen in Figure 7b. Here, a rarefaction from the decreased wall pres¬ 
sure is shown propagating forward to the detonation front. An actual 
overdriven or strong detonation must have an external force driving the 
rear wall at a velocity greater than the CJ product velocity. Fickett 
and Davis [5] give as an example for an overdriven detonation the case 
where another detonation, stronger than the one being studied, is driving 








Pressure 







18 


the rear wall. 

A second possibility would be the case where the wall velocity is 
exactly equal to the product velocity at the CJ condition. Here, the 
detonation pressure is the CJ value and the flow following the front is 
sonic with respect to the front. Figure 8 illustrates case 2. 

The third, and most important case for our study of DDT initiated 
by flame in a closed tube, is when the wail velocity is less than the CJ ve- 
1 city. For this case, the front still propagates at the minimum detona¬ 
tion velocity allowable by the conservation equations, the CJ detonation 
velocity,!^, and the combustion products leave at velocity , sonic with 
respect to the front. Since the wall velocity is less than the product ve¬ 
locity, a rarefaction from the rear of the combustion zone to the wall is 


- 1-- l 


u’evei.upcu. 


This is an iscHtropic expansion characterized by the ratio of 


specific heat, y, assumed for the products (see discussion in previous sec¬ 
tion). Figure 9 illustrates the case where the wall velocity is less than 
the product velocity. As an example of the isentropicity cf the flow be¬ 
hind the front, the pressure ratio between the final reaction state (CJ) 
and the location labeled k in Figure 9 is obtained from the expansion 
wave 



l + O^r ) M 


T2y/(Y-1) 

k J 


where M. is the Mach number of the flow at location k. 
k 


( 21 ) 


1.5 Steady and Unsteady Deflagration 

Along with the three detonation solutions of the Rayleigh line and 
Hugoniot curves (i.e, CJ, strong solution, and weak solution) discussed 
in the previous section, there are ether steady state solutions possible. 
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Figure 10 illustrates a solution to the conservation equations on the lower, 
or deflagration, branch of the Hugoniot curve. For this solution the 
flame propagates through the unburnt reactants at a low Mach number (M«l) 
and as can be seen on the figure, the deflagration is characterized by a 
drop in pressure and rise in specific volume across the reaction zone. Re¬ 
call, for the detonation solution there was a rise in pressure and d„op in 
specific volume. An example of a steady state deflagration would be the 
case of a flame burning on a bunscn burner. Here, the flame is stationary 
and the reactant gas is flowing through the flame sheet at a low velocity. 

When discussing DDT in a packed porous bed, one must clarify exactly 
the fluid dynamics of the deflagration branch. At time t=Q, one end of 
the packed bed is ignited and generating gas from the propellant surface. 

As time progresses, it is the hot gases confined in the region behind the 
ignition front which in a sense drives the ignition front, through the bed. 

That is, unlike che pressure drop across the steady state deflagra¬ 
tion solution, there is a pressure rise across the ignition front. The 
increase in pressure is a result of the gas generation in the region be¬ 
hind the ignition front being confined to a finite volume. Also, in the 
analysis made of the steady state deflagration, all reaction was assumed 
to take place across an infinitely thin re^rcion zone. I* is obvious 
that for the two-phase deflagration problem this is not the case. Be¬ 
cause of the finite zone of reaction, the problem of the deflagration 
flame is at most quasisteady and it may be tha* the nonsteady terms in 
the conservation equations should not be neglected. Because of this 
inherent unsteadiness of the deflagration phase of the DDT problem, the 
reader should not confuse it with the steady state deflagration solut'on 
shown in figure 10. Later it will be shown that if transition from 
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deflagration to detonation does occur, the reaction zone must collapse 
and a detonation solution can be analyzed as being steady and by the 
jump conditions being properly satisfied. 

1.6 DDT in Two-Phase Reactive Flow (Experimental Data) 

As discussed in the previous section, the transition from deflagra¬ 
tion to detonation in a porous reactive medium is an unsteady process. 

Hot gases generated from the propellant surfaces are driven forward into 
the unburnt solid matrix by the pressure gradient developed at the igni¬ 
tion front. This phenomena is not found when a nonporous solid deto¬ 
nates, since only pressure disturbances can be propagated ahead of the 
ignition front. It is this convective heat transfer to the unignited 
propellant and the extended deflagration reaction zone which makes DDT 
in a porous bed unique when comparing it to DDT in either an all gas or 
an all solid regime. 

Much experimental work on DDT in porous material has beer, carried 
out in the past decade. Bernecker and Price [8-10], have published the 
most recent results on DDT in a series of three papers. Other experi¬ 
mental studies prior to these include the work of Griffiths and Grocock 
[11] and Taylor [12]. 

The work presented in Reference 9 by Bernecker and Price is a study 
on DDT in RDX (cyclotrimethylenetrinitramine), a shock sensitive high 
energy explose. In their experiments the RDX was packed into a thick 
walled tube having an inside diameter of 16mm and being approximately 
300mm (12 in.) in length (see Fig. lia [8]). Both ends of the column were 
closed and ionization probes were located throughout tne bed to trace the 
ignition front locus. The RDX was packed in an inert wax mixture and 
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had a mean particle size of 200ym in diameter. Figure lib illustrates the 
DDT mechanism on a distance-time plot for the 91/9 RDX/wax granular charge 
[9]. 4.s shown in the figure, the convective flame front obser/ed in their 

work had no acceleration up to the detonation transition. 

Their experimental work showed this convective ignition front traveled 
at subsonic velocities [0.3-0.9mm/usec) for most cases when the initial 
porosities, <t> o> ranged between values of 0.1 to 0.3. The lower data points 
on the plot represent a post convective compressional wave in the burning 
region which overtakes the ignition front. The length of porous propel¬ 
lant it takes for the compressional wave to overtake the ignition front 
and transit into a detonation is i^j, the run-up length. 


1.7 Topics to be Addressed 

The work nresented in this chanter reviewed those features which dis- 

A A 

tinguish DDT in a porous reactive medium from DDT in other media. Since 
the work that follows in s^usequent chapters will attempt to model this 
phenomena, it was important that the properties of an actual detonation 
be understood(Sections 3-5). 

It was also pointed out that the DDT process being studied cannot re¬ 
present the transition between the lower CJ point and the uppor CJ point. 
That is, because the flame propagates from a closed end and the deflagra¬ 
tion reaction zone is never of infinitesim 1 thickness, one cannot apply 
steady state jump conditions for the deflagration phase. 

The chapter that follows reviews the past DDT modeling efforts and 
inidcates how this work builds for the present effort. Numerical solu¬ 
tions to the flow equations are carried out am solutions are presented 
in Chapter 4. These results will show that the detonation (steady state 
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supersonic wave) is a solution for a certain class of problems. Such re¬ 
sults have yet to be presented by others. 


I 
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CHAPTER TWO 

THE MODEL AND ANALYSIS 


2.1 Introduction 

The analysis that follows attempts to model the situation in which 
a bed of tightly packed granular propellant is ignited at one end. Both 
ends of the packed bed are considered closed, thus, modeling the problem 
discussed in Section 1.4 where the velocity of the rear wall was less 
than the CJ product velocity. Recall for this case that after DDT, the 
detonation propagates through the unignited region at the CJ detonation 
velocity and detonation pressure and is followed by an isentropic expan¬ 
sion of the gases to the pressure at the stationary wail. The gas sur¬ 
rounding the particles at the initial time is considered inert, and to 
be at atmospheric pressure. It is also assumed that the inert gas will 
fully mix with the gases being generated from combustion of the propel¬ 
lant as time progresses. Figure 12 is a schematic of the propellant bed. 

For numerical simplicity the propellant particles are assumed to be 
unisized spheres. Particles of interest range in diameter 50pm < d^ < 200pm. 
T o treat mult -sized particles, one would require _N independent equations 
of mass, momentum and energy for the solid phase, whore N is the number 
of initially diffcrent-uized particles. A solids loading of 74% is the 
tightest possible for uni sized spheres, obtainable by arranging the spheres 
in a face centered cubic. However, assuming granular deformation occurs 
under high stress loads, as discussed by Kuo [13j , greater solids loadings 
may be predicted without error. Obviously, the spherical geometry must 
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he altered for this to occur. 

As the small fraction of propellant particles ignited at time t=0 
burn, hot gases are generated as a function of the pressure-dependent 
burning rate law and surface-to-volume ratio (3/r ) of the spheres. These 
hot gases generated are convected forward through the lattice of unburned 
propellant and flow gradients develop, as dictated by the solution of the 
conservation equations and the necessary constitutive relations. 

Heat transfer from the hot gases to the unignited propellant, par¬ 
ticles, dependent on the velocity of the gas relative to the particles 
and several gas properties (i.e. viscosity, thermal conductivity), trans¬ 
ports energy from the gas to the solid phase. Subsequent ignition of par¬ 
ticles further down the bed is assumed to occur when a critical solid phase 
internal energy is reached [2], This energy can be expressed as a critical 
increase in solid phase temperature, T , since the specific heat of the 
solid is assumed to be known. 

As time progresses the gas pressure behind the ignition front increases 
due to the confinement of the gases from the closed rear boundary and the 
pressure-dependent rate of mass generation in the gas phase. Under certain 
conditions* the pressure gradient can develop into a shock front which 
overtakes the ignition frcnt propagating through the bed. When this occurs 
the ignition front experiences the transition from deflagration to detonation 
discussed in Section 1.6. 

At the transition point the ignited region (zone of gas generation) 
narrows in width and is followed by a region of all gas where the propellant 


These conditions will depend upon the solid chemical energy, granula¬ 
tion (size and loading), burning rate, ignition energy, etc. 
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particles are completely burned out (see Fig. 13). The thickness of the 
reaction zone is a function of the initial particle size and solids load¬ 
ing, and thicknesses approaching 1mm may be possible as r^ + 0 and 4> 0 -*■ 0. 


2. 2 Assumptions 

In order to numerically model DDT in two-phase flow while retaining 
the physics of the problem, several key assumptions had to be made. These 
assumptions are similar to those made by previous investigators (Ref. 2 


and 4 ). 


(1) Both the solid and gas phases are independently treated as con¬ 
tinuums requiring their own conservation relations. 

(2) Each phase interacts with the other. This is modeled by the 
mass, momentum and energy interaction terms in the conserva¬ 
tion equations. 

f3) All propellant particles are unisized spheres. 

[4) Ignition of a propellant particle is obtained when a critical 
energy, expressed as a particle temperature, is transferred to 
the solid. 


1Iant particles are initially surrounded by an inert 


gas at temperature, T 

*o 

(6) During combustion of the propellant, the gaseous products mix 
with the inert gas described in assumption (5). 

(7) Both ends of the bid are closed allowing no gases to escape. 

(8) The specific heats at constant volume, , for both phases are 
constant. 


(.9) When the solid phase, at a given x-location in the bed, displays 
a porosity <p >0.95 and a particle volume one-tenth of its initial 








volume, it is burned out and no longer generating gas. Results 
show this phenomena to proceed smoothly from the left boundary 
and thus not leaving any 'holes' in the continuum. This assump¬ 
tion was necessary to prevent a singularity from arising as r^ -*■ 0 
and T -*• 00 . 

(10) All the product gases obey an assumed nonideal equation 
of state. 

(11) The solid particles are compressible, without heating up, obey¬ 
ing a modified Tait equation of state. 

(12) Cnee ignited, the particles are assumed to burn on the outer 
surface only, at a known pressure-dependent rate law. 

(13) At some initial time, a "narrow” region at one end is ignited, 
burning at the low pressure prescribed. 

2.3 Governing Equations 

Numerical modeling of the two-phase reactive flow process in DDT in¬ 
volves the conservation of mass, momentum and energy per unit volume in 
both gas and solid phases. This is a system of six conservation condi - 

+■ i rtnc uVi i r>Vi 4 - /-\ r*»n n cat n ■£" nnnl i naer Vm fv* av-Uri 1 i r> -I o 1 n 1 

***v**j * vim w -Jvw w*. uv/aaw put umoi uiibiai CVJUO~ 

tions coupled by the interphase mass, momentum and heat transfer terms. The 
conservation equations for two-phase reactive flow have been developed 
previously and References 3 and 4 will provide the reader with the defi¬ 
nitions, assumptions and expressions for the six following field-balanced 
conservation eauations. The e are: 

Gas Continuity 
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Particle Continuity 
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Particle Energy 
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at 


dx 
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(27) 


Here, the relations for the total internal energy in each pnase are 

E <t 5 C v g T g * I “* a " d V V P ‘‘* "p 

The subscripts g and p denote gas and particle ■ espectively. In Equations 
22 to 27, the phase densities, and p^» are defined as 


P ! = p g* 


and 


P 2 = (l-4>)p p 


The porosity, <J>, is defined as the ratio of the instantaneous gas volume 

Hence, the solids fraction is (!-$)• 


to the mixture volume. 
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In addition to the six conservation equations, three constitutive 

relations are needed in order to solve for the nine unknown variables; 

p,p,u,u,T,T,<t>,P and P . These relations include state equa- 
g P g P g P g P 

tions for both gas and solid states and a stress-resistance relation for 
Pp. Appendix B gives a complete listing of the relations used. 

2.4 Improvements 

Since the work reported in Reference 3, certain "improvements" in 
the modeling effort have allowed solutions which may be considered to be 
actual detonations. These improvements are discussed in some detail later 
in the text, but basically include: 

1. Implementation of the necessary gas phase (nonideal) equation 
of state, to insure that at the CJ (Chapman-Jouget) cc. ditions, 
the isentrope provides for a "gamma law" suitable at the hydro- 
dynamic CJ state. (AppendixC presents a review of such an 
equation of state.) 

2. Implementation of a new gas-particle friction coefficient, a: 
developed by Wilcox and Krier [14], for flows at the high Reynolds 
numbers encountered in the developing DDT flows. Previously such 
coefficients were based on data only available fos moderate Reynolds 
number ranges. 

3. Utilization of a modified numerical integration scheme, which 
allowed for a reduction in the grid spacing (and hence reduc¬ 
tion of the time increment) without the usual penalty of exces¬ 
sive computation costs ana numerical instability that often follows 
when the total number of integrations is significantly increased. 
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2.5 Equations of State: Constraint Due to the Detonation State 

As mentioned in Chapter 1, a nonideal equation of state must be uti¬ 
lized for the product gases. The analysis presented here uses a nonideal 
equation of state for hard spheres developed by S. J. Jacobs [15]. Pre¬ 
vious to this study, a covolume-type state equation with data made avail¬ 
able by Cook [16] was used [see discussion in Ref. 2). 

The hard sphere equation of state takes the form 


1.0 + bp + c (bp) + . . . 


(28a) 


where the constants b and c are determined by the value of the gamma law 
coefficient, y, for the product gases as discussed in Appendix B. As 
stated, y is the negative logarithmic slope of the lsentro^e tangent 
to both t’ne Rayleigh and Ilugonrot lines "t. the CJ point m the detona¬ 
tion state. That is, the slope of the isentrope that the product gases 
expand along in the product state. The reader may refer to Appendix C 
for the complete solution of the constants b and c. 

Values for y range from two to three for detonating high density 
explo ives [15]. For the baseline case considered in this study, a 
value of y = 2.05 was selected and the corresponding nonideal equation of 


state is: 


~r~ = 1.0 * 2.5 p - 0.50 p 2 (28b) 

RT g g 

When the above coefficients (2.5 and 0.50) were altered to treat a case 
for y = 3.0, excessively high gas and particle temperatures were predicted 
as one mignt expect. In addition, during the numerical integration severe 
oscillations in ;as and particle temperatures occurred in most cases when 
y > 3. Since one is always constrained by the numerical integration schemes 
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that are employed to handle very severe gradients in the flow, it is under¬ 
standable that previous efforts in DDT modeling [2, 3], which without 
knowing were utilizing high y values for the product gases, almost always 
ran into numerical failure. 

Evaluating the covolume state equation previously used in Ref. 3, an 
approximate value of y = 3.6 is calculated. This extreme value for y may 
be one important reason why the calculations, as reportel in Ref. 2-4, 
were unable to handle the high pressures associated with steady state 
detonations. 

In the solid phase the particles obey a modified Tait equation allow¬ 
ing for compression of the granules. This is written as: 


0 




K 


1/3 


+ 1 


(29) 


* r o *- o -» 

wnere K q is the bulk modulus. Reference 3 discusses the Tait equation in 
more detail. A typical value is K q = 1.38 GPa (2.0 x IQ 5 psi). 


2.6 GJ (Detonation) State 

In the text by Fickett ard Davis [5] equations are presented for es¬ 
timating detonation velocity, D . , and the detonation nressure. p_. . for 
a steady state detonation. These are given for a single phase explosive: 


and 


[ CJ = 2Cr ‘ 1 - )P p 0 e chem 


d cj= e chem 


(30) 


(31) 


Here, p is the initial solid material density and E is the chemical 
p 0 CHEM 

energy liberated by burning the solid. Since the problem being considered 
is two-phase (solid-gas), Eqs. 30 and 31 must be modified to account for 
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this by converting p , the initial solid density, to p~ , the initial 

P 0 o 

solid p hase density 3 P„ (1-40 • Foi an initial solid density p^ = 1994 


o * o 

=0.30, a chemical energy E 
2.05, Eqs. 30 and 31 give respectively: 


kg/m , an initial porosity 4> rt =0.30, a chemical energy E 
and assuming y 


=5.48 MJ/kg, 


and 


CJ 


CJ 


16.06 GPa 


5.92 -HU 
qsec 


These equations were developed from the jump equations for one-phase flow 
where the equation of state for the product gases was assumed ideal. Be¬ 
cause of this, these equations should only be used to get a good esti¬ 
mate of the detonation pressure and velocity. 

In addition to Equations 30 and 31 another method was presented by 
Fickett and Davis [5] for estimates of detonation pressure and detona¬ 
tion velocity. This is hamlet's Short Method and was developed by hamlet 
and Jacobs [17], in the CJ state they are: 


0.5 


P CJ = 

P 


1000 


(32) 


and 


0.25 


°CJ- Atl * Bo o )n "W • 31 - 62 l3i) 

P 

where the constants A = 2.23 (m-kg-s ^(mole MJ) B = 0.0013 

(m 3 /kg) and K = 0.762 (Nm^ kg *(mole MJ) , Again, p is the initial 

3 P 

solid density in kg/m and N the reciprocal of the gaseous molecular 

weight in mole/kg. Equations 32 and 33 give detonation pressure and ve¬ 
locity for a detonating solid explosive going to all gas in the product 

state. To modify Equations 32 and 33 for the two-phase problem,p , the 

D 

solid density, is again converted to the solid phase density P 0 = (l-4> )p 

Pc 
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It should be noted that Equation 33 expresses detonation velocity as 
a function of propellant density, while Equation 31 shows it to be inde¬ 
pendent of density. Like Equations 30 and 31, Equations 32 and 33 should 
only be taken as approximations. 

Using the same input as above. Equations 32 and 33 give respectively 


CJ 


CJ 


» 16.52 GPa 


= 7.83 


mm 

usee 


2.7 Gas Permeability 

One of the key constitutive relations required in the analysis is 
the gas-particle (interphase) viscous force, which governs hot gas pene¬ 
tration into the unignited region of the granular material. As presented 
in Appendix A: 

U 

V = —^ (u - u ) f (34) 

4r 2 g P Pg 

P 

where f is the drag coefficient. Until recently, the packed bed corre¬ 
lations by Ergun or Kuo and Nydegger (as reviewed in Ref. 14) were uti¬ 
lized for f . Thus, the modeling efforts presented in Ref. 3 and 4 used 
the expression of Kuo and Nydegger [18]: 

i2 .81 


f . {276 + 5 (-—-) } 


Pg 




l-d 


(35) 


Equation 35 was developed for 460 < Re < 14,600. Here, Re is the appro¬ 
priate Reynolds number, defined as: 


Re - [(<t>u g )p g 2r p ]/Pg (36) 

Based on experiments at both high gas velocities and high Reynolds numbers, 
Wilcox and Krier [14] developed the correlation: 
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f = 5.06 x 10 S r Re/u 2 
Pg P 8 


(37) 


where u is given in m/sec t r in meters, and the constants 5.06 x 10 
g s P 

must have units of m/s 2 . This expression is not valid for either very 

low gas velocities, since Equation 37 would give f ■* 00 as -*• 0, or 

for very high gas velocities*, since f -*■ 0 as u -*• ®. The equation 

3 5 

was found to be fairly accurate for 10 < Re < (2 x 10 )and 15 m/s< u < 150 m/s. 

g 

While a straight forward comparison is not easy, the difference between the 

value for f as predicted by Equation 35 versus that by Equation 37 can 

4 

be seen in the example. At Re = 10 and $ = 0.4, Equation 35 gives 

f = 53,600 and Equation 37 gives a value of f = 5448. To utilize 
Pg H Pg 

Equation 37 one must specify the average gas velocity, u , and the par- 

g 

tide radius, r o , that were used to obtain the Reynolds number. A kine- 

“6 2 

matic viscosity y /p = 1.8 x 10 m /sec, a particle radius r = 1.0mm, 
g g P 

and an average gas velocity u^ = 30.5 m/sec were used to obtain a Reynolds 
~ 4 

number Re = 10 . From this particular example, the Wilcox/Krier correla¬ 
tion (Eq. 37) allows about 10 times the permeability, i.e., 1/10 the viscous 
drag force as correlated by the Kuo/Nydegger relation. 

9 'll Hi'S nnc tn flow nr-nr <»<;«; 1 uadi no fn DOT di rn<; «tA<' in Tnllnw- 

-- — --- r - ---© - * ---- - -- - 

ing chapter, clearly indicate that sufficient gas permeability is necessary 

to allow for a detonation transition. 


•Data from Ref. 14 was limited to u < 


300 m/sec. 
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CHAPTER TmEE 

NUMERICAL INTEGRATION 


3.1 Finite Difference Mesh 

To solve the Eulerian formulated system of conservation equations 
discussed in Chapter ? with the constitutive relations listed in Appendix 
B, the length of the bed being integrated over is divided into I segments, 
each a constant Ax in width (i.e., x. = jAx; j = 1, 2, 3.,.. I). The value 
required for Ax will be discussed later in the chapter. 

At t = 0 values of the nine independent variables; p , p , u , u , 

* 8 P 8 P 

T , T , P , P and $ are initialized at each ith x-location in the erid. 
g P g P 

Before incrementing the primary variables (i.e. mass, momentum and energy) 
to the future time, t - t Q + At, the auxiliary variables in the equations 
(e.g.,drag, gas generation, heat transfer) must be computed at the present 
time, t = t . The nine equations are then solved at the incremental time, 
t = t + At,by a modified Lax-Wendroff finite differencing scheme. This 
method along v .th another tested are presented in Appendix D. The second 
method, developed by Rubin and Berstein [19], was implemented for several 
test cases and proved to give results similar to that of the Lax-Wendroff 
scheme. However, to keep all results consistent, the Lax-Wendroff scheme 
was cho’-en ?.s the integrator. 

The time Increment, At, over which the equations are solved is cal¬ 
culated by the Cour&nt, Fredrich;;, Levy stability criteria [20], for hy¬ 
perbolic equations, 

XAx 

s (Fi^n 


At 


( 38 ) 





1 
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In Equation 38, the term c is the mixture sound speed and |u| is the maxi¬ 
mum gas velocity in the bed. Also. X is a stability constant, less than 
unity. For most cases X = 0.5 was used. This smaller time increment 
was necessary in order to integrate the equations ’.'hen large gradients 
developed in the flow. 

3.2 Initial end B oundary Conditions 

To initialize the problem, the bed is assumed to be quiesent, i.e., 
at a constant gas temperature and constant gas pressure throughout the 
length of the bed. The spherical propellant particles are typically fix¬ 
ed at a constant solids loadings, although a variable initial porosity 
can be treated. Then to initiate the flow, the propellant at the first 
few grid points is assumed ignited and generating gas. To be consistent 
with the fiat initial pressure profile, all gas and particle velocities 
at time t = 0 are set equal to zero. 

In the past an exponential pressure profile was constructed at the 
initial time in order to speed up the deflagration process. This was 
inaccurate, however, in that it was unsupported (i.e. , no fluid motion 
was associated with the pressure gradient). Nevertheless, the computa¬ 
tion time necessary for the pr 'ssure gradient to develop on its own was 
too great and the initial pressure gradient was implemented to speed up 
the proces-' Since then the computer code was modified to include an in¬ 
tegration process which cuts the computation .ime drastically and allows 
for the use of a more realistic flat initial pressure profile (P ~ P ). 
Table 1 is a summary of typical input data for the cases studied. 

3. 3 M odificat io n to the 1 n t egrati on Scheme 

To solve on a digital nputer the flow which is to represent a DDT 
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Table 1 

TYPICAL INPUT DATA 


__P ARAMETER _ 

Burning Rate Index 

Burning Rate Proportionality Constant 

Initial Bed Porosity 
Particle Diameter 
Bed Length 
Grid Spacing 
Chemical Energy 

Gas Specific Heat 


VALUE 


0.8 < n < 1.2 

b = 0.001 — [~1 

sec l PSI J 


0.25 < 4> < 0.50 

— o — 

50um dp <_ SOOum 
L = 25.0 cm 


Ax = 1.27 mm 

4 — < If < 7 Md 

kg CHEM 7 kg 

If T 

L 0 kA^ c v ^ L9 


KJ 

kg^k 


Initial Bed Temperature 
Ignition Temperature (Bulk) 

Ignition Energy 


T 294 °K 

g 

T. = 303 °K 
ign 

E. => 

ign 


9.0 


KJ 

H 
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phenomena, requires the repetitive integrations of many equations. All 
of the conservation equations and constitutive relations are solved indi¬ 
vidually at each grid point i\i the entire bed, for each time increment. 

A study of several test cases showed, for instance, that when the ig¬ 
nition front was at a given x-location, there was little activity several 
grid points ahead of it. That is, the gas and particles were at velocities 
close to zero and all transport coefficients were negligable in the region 
not far in front of the ignition front. This phenomena occurred for both 
the detonation state where the front was propagating between 5-8mm/usec, 
and the deflagration state where the propagation velocity was much less. 
Because of this, an addition was made to the computer code which allows 
the code to only integrate the active region and bypass the inactive zone. 

Ahead of the ignition zone, the computer code locates the nearest 
point to the zone where there is no significant particle or gas movement. 
This point is then designated as the new front boundary of the integration 
region for that particular integration step . When the pressure front 
builds the ignition front moves rather rapidly through the bed (5-10-™^.-) 
and a new integration boundary must be located after each time increment. 

Figures 14 and 15 show s comparison of a case run without the moving 
integration zone (Fig. 14) and one with (Fig. 15). To assure a correct 
solution the actual front integration boundary is extended a few grid 
points beyend the location calculated by the code. The addition of this 
new logic to the current code reduced computation ti e by at least one 
half, for the same time increment and grid space. 

3.4 Artific i al Smoothing 

Inherent to the solution of the system of interdependent conservation 















equations is a numerical instability. A small perturbation can in some 
cases amplify with each time increment and eventually destroy the numeri¬ 
cal solution. This phenomena can start at the first time increment and 
in ten to twenty integrations the oscillations can be so large that the 
solution becomes unstable and terminates. In order to smooth out these 
oscillations before they amplify, an artificial smoothing routine must 
be incorporated into the code. 

From experience in integrating the two-phase flow equations (Eqs. 
22-27) with significant, nonlinear source-sink terms, the problem of nu¬ 
merical instability occurs often enough to warrant artifical smoothing. 

An extensive study by the author and other investigators [2-4], shows the 
final solution to be independent of any artificial smoothing used- 

The analyzation of numerous test cases has shown, suprisingly. that 
the stable solution did not require smoothing of all the variables. It 
is obvious that the following results have an inherent dependency on 

smoothing . Extreme care has been taken to minimize these effects on the 
qualitative trends and quantitative results predicted. Nevertheless, it 
should be obvious that smoothing techniques can supply variability in the 
predicted parameters which are reflections of the scheme utilized, and not 


necessarily of the conservation equations. 

The particular artificial smoothing technique found to be useful is 
a three point averaging technique where the variable v determines what 







4? 


whore 



For example, if v = 0.1 then W retains 80% of its original value and is 
influenced by 10% of the U values on each side. This is shown graphically 
in Figure 16. 

For all cases studied, a minimum value of V was used that would pro¬ 
vide for a stable solution. This required repeating each test case nume¬ 
rous times, lowering v ea?h time until the solution went unstable. 

3.5 Grid Spacing 

It is obvious that in order to minimize computation time, a maximum 
value of Ax which gives a stable solution should be used. Based on the 
calculations carried out in References 2-4 a grid space of Ax =■ 1.27 mm 
was shown to be the largest Ax that would provide a solution to a flow 
problem which exhibited rather steep gradients. 
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CHAPTER FOUR 


RESULTS COMPUTED 


4.1 I ntroduction 

This chapter will present the calculations made on the possibility 
of DDT occurring in packed beds of high energy, granulated, unisized pro¬ 
pellants or explosives. It is obvious that there are an infinite number 
of loading combinations possible for the DDT study. Fortunately, the work 
of Hoffman and Krier [3] ana Krier and Gokhale [2], in which conditions 
of rapid convective flame spreading have been calculated, is available 
and can be used as a starting point. It should be noted that none of 
the calculations made in References 2 or 3 predicted DDT. As pointed out 
in the previous chapters, there are probable reasons why this has not been 
accomplished and it is expected that certain improvements and modifica¬ 
tions will now allow for the calculation of the steady state detonation 
solution. 

The study made of a DDT potential attempts to model conditions simi¬ 
lar to the test conditions of the Bernecker and Price work [9] fi.e., a 
long column of granulated material in a closed pipe ignited by an ener¬ 
getic ignition material at one end]. In order to model the flow tran¬ 
sient, one must assure that the lengtn of the bed exceeds £ , the run-up 

length to detonation. Since the experimental work of Bernecker and Price 
[9] indicated that a 10 in (25 cm] bed was sufficient for most of their 
experiments where DDT occurred, this length was selected as the longest 
bed length to be considered. 
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Since unisized spheres are being treated, the initial porosity can 
be no less than $ = 0.26, although randomly packed unisized spheres 

generally give a high porosity, about = 0.40. Therefore, for this 
study 0.26 <_ < 0.40. The initial particle radii studied were also in 

the same range as those considered by References 2 and 3. Results in 
these studies showed that particles must be less than one millimeter in 
diameter in order to generate sufficient gases for the rapid flame spread¬ 
ing phenomena. 

The chemical energy of the material considered is in the range of 
explosives or high energy propellants of the nitramine family (i.e., HMX, 

RDX). Thus, the chemical energies studied were always larger than 1000 
cal/g (4.15 MJ/kg). Other parameters one must consider in the DDT studies 
are the burning rate properties of the propellant. Again, the values used 
in References 2 and 3, which attempted to model the burning rate of an 
HMX solid propellant, were utilized. However, the burning rate index, 
n, is a parameter which is explicitly studied. 

In this analysis, the deflagration will be initiated by assuming at 
time, t = 0, that except for a small portion of the bed at one closed end, 
the bed is quiescent and unreacting. As has been mentioned in Chapter 3, 
a closed end situation is considered and hence at the two end points ( x = 0, 
x = L) it is assumed that all flow gradients are zero and that the veloci¬ 
ties of the particles and gas must also be zero. 

The following section begins with a solution that indicates, for the 
first time, a steady state detonation can be predicted. This result is 
then compared to conditions where no such detonation solution occurs. 
Additional calculations will indicate the sensitivity of the initial 
porosity, <J> , chemical energy, E , burning rate index, n, and ignition 

O wiitM 
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energy. 



on the run-up length to detonation. 


4.2 Calculations 

Figures 17a and 17b present a case where a transition from deflagra¬ 
tion to detonation has occurred. For this example, the burning rate in¬ 
dex was n - 1.0, the chemical energy E^^ = 5.48 MJ/kg and the particle 

radius r = 127 um (.005 in). 

P 

The pressure-distance profiles at five separate times after ignition 
of the propellant at x = 0 are shown in Figure 17a. Examining the pro¬ 
file for t = 50 ysec, one can observe that the profile is characterized 
by a shock front at x = 23 cm followed by a smooth expansion back to the 
wall at x = 0 cm. The pressure in front of the shock is at atmospheric 
conditions and, therefore, negligible with respect to the pressure behind 
the shock. 

The ignition locus plot for this particular case is shown in Figure 
17b. Here, the ignition front moves through the bed at. a low subsonic 


velocity for the initial ten microseconds and then accelerates to reach 
a steady state velocity of =7.2 mm/ysec. This occurs within 12 to 15 

cm from the ionitprl f»nd . Note that th^sp df*tnnat i on solutions for P 

and are in fair agreement with the approximations made in Ch? 2 
(Eqs. 30 and 31) for the given input » V Yj)- 

Obviously, hydrodynamic steady state analysis (like Eqs. 30 and 31) 
cannot guarantee that a transition from deflagration to detonation will 


occur. However, it seems that "critical" values of porosity (related to 
gas confinement), gas generation rates, and chemical energy will provide 
for a DDT. For example, when the burning rate index was lowered to 
n = 0.8, as shown in Figures 18a and 18b, the steep pressure front associated 
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with a detonation did not develop. Correspond!ngly ? no detonation speed 
was predicted. In this example the peak pressure in the bed never ex¬ 
ceeded 5 GPa, no shock was predicted, and only a steady convection-dr’ven 
front of 2.2 mm/usec occurred at 100 ysec after ignition of \ = 0. It is 
not clear at this time whether or not this can be defined as a "low-ve¬ 


locity'' detonation. 

According to Equation 30 in Chapter 2, for the case where a transi¬ 
tion to detonation actually occurs, the steady state detonation pressure, 

P , should increase linearly with the chemical energy, L . Also, Equa- 

HIlM 

tion 31 states that the detonation velocity, D^. , is a function of the 
square-root of the chemical energy. Figures 19a and 19b present the re¬ 
sults of a calculation where all parameters were identical to those used 
to give the DDT results of Figures 17a and 17b, except * 6.85 MJ/kg, 

an increase of 25%. The steady state shock pressure predicted for this 
case, shown in Figure 19a, was 21 GPa. This represents a (21/16.4 = 1.28) 
28% increase in pressure over the first case. The predicted detonation 
speed (Fig. 19b), calculated from the slope of the x-t diagram, was 8.70 


nun/usec. According to Equation 31, the ratio of for the case given in 
Figure 13b to that presented in Figure 17b should be 


A6.8S75.48) - 1.32 

This is approximately the increase predicted. 

Table 2 summarizes the detonation pressure, P , and the detonation 
velocity, D. T , (which were the end results of the TUT calculations) all 
as a function of the propellant chemica 1 energy. These predicted condi¬ 
tions are compared with the approximate steady state hydrodynamic solu¬ 
tions discussed earlier, i.e., Eqs. 30 and 31. The excellent agreement 







'Location (cm) % Press 



.g, 18a. Pressure History During Deflagration with no Transition 
(n = 0.8) . 



Fig, 18b. Ignition Front Locus with no Detonation [n = 0.3). 










(cm} 



Fig. 19a. Pressure History Leading to Detonation with 25% Increase 
THEM 


’i: u-onijiiie lo case shown Vii i ig. 17a) 
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of tho detonation pressure with the analytic solution should be noted. 

However, the predicted value for the detonation velocity, D^, from 
the hydrodynamic solution (Eq. 31) shows to be slightly less than the 
value predicted by the code for all cases. Although one cannot judge 
which of the two values, detonation velocity or detonation pressure, is 
more accurate, the percent increase in the detonation velocity as the 
chemical energy is increased compares favorably with the hydrodynamic 
solution. 

4.3 DDT Run-Up Length 

The run-up length to detonation is defined in this report to be the 
distance from the closed end where the bed is ignited to the location 
where both the peak pressure and the detonation speed are constant, i.e. 
the equilibrium steady state solution. 

Figures 20a and 20b plot the predicted run-up 1 
as a function of the burning rate pressure-index, n (Fig. 20a), and the 
initial bed porosity, 4> Q (Figure 20b). The "no-solution" boundary indi¬ 
cates that, with the integration scheme used, the mesh size would have 
to be drastically reduced (thereby decreasing the integration time incre 
ment) in order to obtain a stable solution. This is a costly exercise, 
but future work is planned to increase the solution regions. Figure 20a 
clearly indicates that, as expected, one cannot achieve detonation if 
che burning rate during the deflagration phase is not sufficiently large 
(see "no transition" boundary). For the solids loading considered 
(1- <J) 0 - 0.70), it would appear that a minimum DDT run up length is 5 cm 
for particles of 250 pm in diameter. 

Figure 20b begins to resemble the required "U-shaped" curve of l ( . , 
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TABLE 2 

Comparison of the Predicted Detonation State with 
an Approximate Hydrodynamic Solution 


e chem 

V ’ CJ (predict 1) 

P CJ (Eq. 30)* 

D (predicted) 

L J 

D CJ (Eq. 31)** 

£ CJ 

4.11 MJ/kg 

14.0 GPa 

12.0 GPa 

- 

6.48 -22- 
psec 

r , _ mm 

5. Iz - 

psec 

89 mm 

5.4S 

16.6 

16.0/ 

7.25 

5.92 

51 

6.85 

21.8 

20.08 

8.24 

6.62 

44 

8.22 

24.6 

warn 

^ - 17 | / . 2 j 

i 

7 O 

•J o 

10.96 

30.0 

33.2 

10.58 ! 3.37 

! 

■ 


n 

u 
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versus initial porosity. DDT experiments by Korotkov et. al. [21] show 
a similar behavior (see Fig. 21). One would expect that for a relatively 
porous bed, 4 > q > 0.60, no transition will occur since local pressure con¬ 
finement is limited. If the porosity is too small no gas penetration 
for the accelerating deflagration wave will occur. The net result is a 
porosity where a minimum run-up length occurs. 

Figure 22 presents a study on the effect of the chemical energy, E rucu , 

LilLM 

on the run-up length, , the values which were presented in 'Iable 2. 

As expected, the run-up length to detonation increases as the amount of 
chemical ener Q y decreases. Again, as expected there is a minimum value 
where no transition occurs, = 3.0 MJ/kg. For this case it does ap¬ 

pear that a constant but small run-up distance is still required as the 
chemical energy increases beyond the values studied here. 

4.4 Detonation Reaction Zone 

Gne measure of the fact that detonation occurs is a plot of the 
reaction zone width versus the locus of the ignition front. This is pre¬ 
sented in Figure 23, where the reaction zone is defined as the region 
where particles are ignited and generating gas (i.e., are not burned out). 
As shown, the zone initially increases during the deflagration phase as 
the convective heat transfer provides energy to ignite more and more of 
the bed. The zone then collapses to a thin (constant) width as the sur¬ 
rounding high gas presssure causes the part'cles to burn out rapidly. A 
steady react!on zone thickness of approximately y mm is predicted. How¬ 
ever, most of the gas is generated in a small region immediately behind 
the ignition front where the particles are still relatively large, .iote 
that in all the cases reported here, the initial particle diameter was 
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Fig. 20b. Variation in Initial Bed Porosity, <p } on Run-Up Length, 
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2.S0 pm. Obviously, smaller particles will provide for a thinner detona¬ 
tion reaction zone. Figure 24 presents a porosity-distance profile for 
a case where DDT occurred. Here, a porosity of <J> = Q.95 is for all in¬ 
tents and purposes the condition when the propellant is burned out, i.e., 
no generation of hot gases. 

4.5 Comments and Interpretations 

The results shown in Figures 17-24 have clearly indicated that, as 
expected, high solids loading of relatively small particle size ener¬ 
getic propellant with perfect confinement will transit into a detonation 
in space domains of several centimeters. One would expect that propel¬ 
lant properties and packing configurations have limits where detonation 
cannot occur and this was clearly shown in some of the figures which show 
the run-up length versus property parameters. In conclusion, it will be 
useful to review these studies to determine the properties and configura¬ 
tions which minimize a DDT hazard. 

For example, for a fixed chemical energy, E^^ = 5.43 MJ/kg, a fixed 

ignition energy, Ah. = 9.0 KJ/kg, a fixed particle radius, r ■ 127 pm 

^ k ^ 

tO.0,05 ini and a fixed solids loadinv. -fl-rtil - 0.7. the burn!nv rate index. 

, must be larger titan n =* G.St> if DDT is to occur. This was shown in Figure 
20a. Of course, had the burning rate coefficient, b, been a different num¬ 
ber, this exponent may have been different. A general statement >.n then 
be made that the burning rate, and hence rate of gas generation, should be 
kept as low as possible to minimize DDT. Conversely, the higher the burn¬ 
ing rate, the better the chance of a transition to detonation occurring. 

Figure 20b showed the run-up length to detonation versus initial 
porosity, <(, . As was stated in the previous section, for all parameters 




. 24. Porosity-Distance Profile for a Case where DDT Occurred. 
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equal, there is a maximum initial porosity where no transition occurs. 

Since a randomly packed bed of unisized particles has an initial porosity 
of the order <J> o = 0.4, this figure indicates that this is in the region 
where DDT potential is at a maximum. As one reduces the initial porosity 
(that is increases the solids loading) the run-up length to detonation 
slightly decreases until, as experimental work of Bernecker and Price [9] 
has indicated, there is a mimimum initiax porosity where DDT cannot occur. 
Since at these initial solids loadings multi-sized particles and mechanical 
packing are required, these lo lings are not of interest, whereas the ran¬ 
domly packed loadings are of interest. 

A similar comparison of the run-up length to detonation was shown in 
Fig. 22 where the chemical energy content was varied. Recall that a chemi¬ 
cal enerflv. = 4. IB MJ/ke 11000 cal/el. can be considered an ener- 

getic propellant material. Based on the results shown in Fig. 22, one 
can state that less energetic material than this has little chance of en¬ 
countering a DDT. Doubling the chemical energy from E^^ = 4.18 MJ/kg to 
E CHF.M = 8-36 MJAg, represents a reduction in the run-up length of only 
about one half. To summarize, high energy propellant of the nitramir.e 
family where 4.18 MJ/kg, definitely fail within the regime of a 

DDT hazard if properly confined. 

The final comparison of this type is shown in Figure 25. This shows 
the z'un-up length to detonation versus the particle radius. As one would 
expect, there is a maximum particle radius (i.e., surface-to-volume ratio) 
where transition ls not predicted to occur. This is indicated by the "no- 
transition" boundary on that figure. The figure also indicates that as 
the particles get smaller in size, the run-up length to detonation also 
decreases, as expected. Reducing the initial particle size, r , to 

^ n 


Numerical Instabilities 
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Fig. 25. Run-Up Length to Detonation Versus Initial Particle 
Radius, r 
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values less than 25 pm results in gas generation rates, per unit-volume, 

that are so large finer grid spacing must be utilized to assure stability. 

This expensive task has been delayed and is recommended only after improved 

numerical integration schemes have been developed. 

A final topic studied dealt with the effect of ignition temperature 

(or more appropriately AE. ) on the run-up length to detonation. For 

the results shown in this chapter a nominal value of T. = 545°R was used. 

ign 

For this study = 0.03, r pQ = 127 pm and n = 1.0. 

Calculations were made in which T. varied over the range 530°U < 

ign 

T. < 575°R. For a constant initial bed temperature of T = T = 530°R, 
ign - p g 0 r 0 

this represents a range of ignition energy of 0.0 < <_ 25>, 

For most of the ignition temperatures tested there was little change in 


the steady state detonation pressure or velocity. Only when 
approached the improbable value of AE^ = 0-0 KJ/kg did the values 
change significantly. 
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APPENDIX A 

JUMP CONDITIONS FOR TWO-PHASE REACTIVE FLOW 


In Section 1.3 (Jump Conditions) of the text the reader was provided 
with a review of the jump conditions across a shock discontinuity for one¬ 
dimensional, one-phase flow with heat addition. This Appendix will out¬ 
line the development of the jump conditions for one-dimensional, two-phase 
flow with heat addition. As d scussed in Section 2.3, the heat addition 
for two-phase flow comes as a release of chemical energy from the ignited 
propellant particles to the surrounding ga..c;. For a detonation this is 
assumed to occur in an infinitesimally thin reaction zone. 

In order to evaluate the jump conditions for twe-phase flow across 
a steady state combustion wave, the conservation equations (Eqs. 22-27) 
are first written as conservation equations for the mixture. 

Mixture Continuity: 

—- [<f>p u ] + 4 — [(l-$)p u ] = 0 (A. 1) 

ax L g g J dx p p J J 

Mixture Momentum: 


[$0 u 2 + 4P j + T - [(!-<)>)P + (l-4>) P ] 
dx g g J dx P P p J 


= 0 


(A. 2) 


Mixture Ej ergy: 


4- [<J>P U E + <J>U P ] + ~r~ [ (1 -d>) p U E 
dx L ^g g g T " g S dx 11 p Pl 


+ (l-4>Iu P ] = 0 
P P 


(A. 3) 


As in Secti,.-> 1.3. the unreacted or "cold" end is denoted by subscript A 
and the products or "hot" end is denoted by subscript B in the analysis 
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that follows. 

The mixture conservation equations (Eqs. A.1-A.3) are now integrated 
froi the "cold" end to the "hot" end of the bed with the boundary conditions 
<KB) = 1.0 

V A) " V* 1 ■ V A 

V“> • V B 

P S (B) = p b 

Integrated Mixture Continuity: 


P B V B - Vg A V A + C 1 ^A )p p A V A 


Integrated Mixture Momentum: 


p » - * p „ J 

-A 


(A. 4) 


* " V">p a v a ♦ V 


Intcg’.r • tMixture Lb.eigy: 

g T 


, pv vl * V« “ K v v! 


(A. 5) 


(l-itip v(L : v- ; i * 
h 'j j A 


i’jVj " [Cl-tn’ p V 


(A. 6) 


At first glance. Equation A.4 -'o h.. f ;eem to : -.asily evaluated know¬ 
ing the conditions at both ends of the flow. However, as was pointed out 
by Kuo and Suinmerf ield [23], they are not simple ilgcbTaic equations un¬ 
til P , the stress transmitted through the parti-;'- phase at location A. 

F a 

is evaluated. This was stiown to be obtainable by • .tegratlng the momentum 
equation from tin. cold end. A, to the location wi >». e the particles are no 
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longer in contact, denoted by subscript. C. From this you obtain 

fC C 

-P (1-$.) - - V dx + ” u dx 

P A A P 


tC rC 

V dx + ” u dx 

P 

A A 


d u* (1-40 

V 2 (l-<j>) 

Lp p J c L 

p -1 a 


It should be noted that if the flow situation being modeled consists 

of an infinttly long bed, then the particle stress ac location A, P , 

P A 

will also be equal to zero. This is true since pressure disturbances ini¬ 
tiated at the combustion zone (i.e., compression of the propellant par¬ 
ticles) are propagated through the medium at the speed of sound of the 
solid, and would therefore take an infinite amount of time to reach loca¬ 
tion A. 

Since all experimental work consists of finite length beds (the nu¬ 
merical model studied in this report was L = 25 cm)> a comment on the 
integration limits of Equation A.7 is appropriate. If , te region where 
particle drag is of significance extends from the combustion wave (C) to 
the wall condition (A), then the integrals in Equation A.7 must be evalu¬ 
ated after each time increment and hence the particle stress at the wall, 

P„ , is changing. However, if the region of significant particle drag 
^A 

is assumed to be only that z >ne immediately ahead of the combustion zone, 

then the integrals in Eq. A.7 are always constant and, therefore, P is 

P A 

also constant. 

A final comment on the particle stress is that if the time to detona¬ 
tion is less than the time it takes a pressure disturbance to travel 
through the solid matrix to the rear wall, then the rear wall has not 
felt the pressure disturbance and hence, P - 0. 


0 . 
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APPENDIX B 


CONSTITUTIVE RELATIONS 


As discussed in Chapter 2, a system of nine independent equations is 

necessary in order to solve for the nine unknown variables describing two- 

phase flow: p,p,u,u,E,E,P,P and 4>. 

F g P g P g P g P 

Cf the three additional equations necessary for the solution, one is 
the nonideal equation of state for the gas phase as described in the text. 


P 

_JL_ 

RT 

g 


1 + bp 

g 


c (bp g ) 2 + . . . 


(B.l) 


Values for the constants b and c are discussed in Appendix C. 

The second represents an equation of state which relates the solid 
density to the stress on the particle. (See Ref. 3). This is also men¬ 
tioned in the text. 


p 3 K 

3 - ((-2-) - i) ~ 

p “p J J 3 


(B.2) 


Finally, the third additional constraint is a relation for the par¬ 
ticle phase stress, P , as a function of the solids l.-ading and the material 

bulk modulus, K (also see Ref. 3J. 

o J 

As discussed in References 2 and 3, one must also specif> functional 
relations for the following: 


a) 1 ’ = mass generation rate per unit ’'olume 
T = —• (l-it) p r 


(B.3J 


Here, r 

p 


the instantaneous panicle radius and 


is the surface burning 
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rate specified as a function of pressure (and possible partible temperature), 
For all cases run in this study. 


r = b" 


CE.4) 


where r has units of (in/sec), P (psi) and _b is of the order (1 x 30'') 
in/sec . 

(psi) 

b) V = interphase (gas-particle) viscous force (as discussed ir. the 
text). 


where 


V = (u - u ) f 

4r 2 8 P P8 

P 


f = 5.C6 x 10 5 r Re/u Z (see Eq. 37) 
P2 P g 


(3.5) 


c) Q = interphase heat transfer rate 


Q = — (1-40 h (T - T ) 
r p PR g P 

In the analysis carried out here, the heat transfer coefficient was 


(B.6) 


h = 0.65 [Re) U,/ (Pr) U 

P P 


(B. 7) 


where k^ is the thermal conductivity, Re is the Reynolds number and Pr is 


the Prandtl number. 
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APPENDIX j 


NONIDEAL EQUATIONS OF STATE 


As described bv Jac obs | Uj, the nonideal equation of state is ini¬ 
tially written ns a polyru.iniai expansion in gas density. 


Pv 

RT 


1 +• x 


a 

ex" -*• d> 


CC.l) 


Here, x = bp - —. 

v 

If C v> the specific heat at a constant volume, is assumed constant 
over the varying temperatures, as was in the model, Eq. C.l can be 
written as: 


P v = (^-1) E f(v) (C.2, 

C 

In Eq. C.2, V. is the specific heat ratio for the ideal gas limit, y. -- , 

1 v 

and f(v) is the gas density polynomial expansion. 


,., ■ , b ,b. “ 

i (v i - 1 + — + c (—) *• 

' v v 


(C 


For simplicity t-ie higher order terms are dropped. 

The local speed of sound, a, in the detonation state is defined as: 


2 

a 





Here, y is the e ffe ctive specific J.sat ra-'o in tlie detonation s tate. 
For detonating ••• j lusise:- ■ fvp 1 . ■ ’ 1\ we<m > v.. to • h ee . The 

equality in Equation . 4 ■ ■■ no-.- s< t : r< ,<■ 
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Y = £ (—) = £ (IB.) (l£) - i i (IB.) r. .L_) 

Y P C 3p J P '■3v' ’■Sp' 1 P l 3v J *■ '> 1 


s p 


o: 


“V ,3P, 

v = P ( 9v j 


(C.5) 


Since P -- P (fc, v) , 


d! * a!. : - (~) dv 


(C.6) 


and 


dP J»j_, di 

^d v 3E' d\ 


i.\ ,dv, 

'dv _ ‘-dv 


(C.7) 


Substituting Eq . C.‘ : an-- tne constant entropy thermodynamic relation 


P = ■ 


p-.o F.q 


- j 

:T 


' : * P ) 


(c.») 


f: '' -t, t: - ; 


• T Li ..r- * 

I of eq S ? a. 
•7 ) ■ 


M/af •••**• a rm<»ar i n a i ri Fn C - ^ 


v ..s taken holding t constant, 
f’(v) 


•3v 


(C. 10) 


and ‘, conu «rt a *. i. 


r' o t: is taken. adding v 


l-L 


f * 

iM —. — • 


(c;. ] J. : 


m. • si! 


i*. n j j* 


, i> r.». •»* .a> 
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the form: 


y 


(Yi-13 f(v) 


1Y.-1) E f(v) 


CY.-l) E f' (v) v 
Pv 


But, since from Eq. C.2 


CC.12) 


Pv 1 
(Yj-DE " f(v) 

Eq. C.12 is now: 

Y = (Y-l) f(v) - 1 - v f-^p- (C. 13j 

Knowing approximately what value of y an explosive exhibits, values of 
the constant b and c can be fit to equate Eq. C.13. 
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APPENDIX D 


FINITE DIFFERENCING TECHNIQUES 


This section will give the reader a more detailed description of the 
finite differencing scheme referred to in Chapter 3, the Lax-Wendroff 
center differencing technique. 

In order to solve the system of six hyperbolic partial differential 
equations describing two-phase reactive flow, the Lax-Wendroff scheme pre¬ 
sented in Reference 19 for one-phase flow was modified to model two-phase 
flow. This modification was necessary because of the interactive source- 
sink terms involved in the governing equations. 

The six conservation equations described in Chapter 2 (Eqs. 22 to 21) 
can be written in vector form as: 


30 t 3F 
9t * dx 


(D.l) 


In this notation U, F and C are the vectors 


P 1 



P 2 


° 2 U p 

P 1 U H 


(p i% * * y 

P,u 

2 P 

F = 

tP 2 U p * Cl~<}>)F p ) 

P l E g T 


(p.u E + u PI 

1 n T g g 

P 2 E g T 


P 2 U P L P t ♦ (1-<J>) Up P ) 
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C = 


r 

-r 


r u - v 

p 

-r u + v 

p 


r ‘ E »EM * -$> ' ® S ' « 


r ^ E CHEM ^ ^ 13 u p *■ Q 


The Lax-Wendroff solution, as described by Richtmeyer and Morton [22] 
begins as a Taylor's series expansion accurate to the second-order. 

\ 2. t 


,.t+At „t A4 . [Sul* 1 . 2 [3 u 

U = U + At -s— + -=■ At —rr 

x X at I 2 U 2 

*- -i x -ot - 


(I). 2 ) 


In order to make use of Equation V.2, the time derivatives must be replaced 
by derivatives in space. 

3F i 

By using Equation D.l and introducing the matrix A, where A . = ~ 


3 II 3_ 
2 = 3t 

at 


1-1 

t 

X |*TJ 

3C 


’3Fl 

3 C 

3__ 

r 

3U ' 

3C 

3_ [ 3F 

= 3t 

3x 

_3tj 

3 3c " 

3x 

l a 

at _ 

at ‘ 

3x 3x 

. m 


CD.3) 


Replacing x-derivatives with finite difference quotients (i.e., 

3F/3x = (F x+Ax - Ax 3/2Ax), the Lax-Wendroff solution can now be written 


,t»-At 


u* - - A -i 

x 2Ax 


F l . - F* . j + AfC 

x+Ax x-Axj 


[At] 2 

' t 

r t ,.t 1 

. t 

r t t V 

[AxJ 

x+Ax 

L 2 

F . - F 

x+Ax xj 

- A 

| x-Ax 

2 

F - F . 
x x-Ax 


+j At*C (D.4) 


r 
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A*" . denotes A [x . +4 ll*] ■ For A a constant matrix the Lax-Wendroff 

x+Ax 2 x+Ax 2 x J 

£. 

solution takes the form 



CD.5) 


Richtineyer [22] developed a variation of the Lax-Wendroff solution al¬ 
so of second-order accuracy. This method involves two steps and eliminates 
the use of the A matrix found in the Taylor's series expansion solution. 

Richtmeyer's first step involves computing intermediate values in 
time and x-location of U. 


„ At 
t+ -s- 

J * 

X+ Ax 

2 


■ 1 K'x.ix * < ] * m [ r x.Ax - F x ] ' r • c x 



The first step (.Eq. U.t>j is of first-order accuracy. However, the final 
step (Eq. D.7) is accurate to the second-order since r quantities of 
8 (Ax) are differenced over Ax. For the linear case where F(U) = AU 
Equations D.6 and D.7 combine to yield the Lax-Wendroff solution, Equa¬ 
tion D.5. 

A two-dimensional (x-t) grid is helpful in understanding how the 
Lax-Wendroff differencing schome can be applied to the numerical solution 
■>\ the conservation equations on a digital computer. In Figure D.1 a known 






S2 


value of the vector U (mass, momentum and energy per unit volume) at a 

time t = t and located at x * x q is indicated on the grid by a solid 

dot at (x q , t Q ). Locations where the vector U is yet to be solved at 

are indicated by open circles. Figures D.2 and D.3 illustrate the 

Richtmeyer two-step routine using the x-t grid. For instance, in Figure 

D.2 the vector U at (x - . t * 4^-1 is solved for by using known values 

O 2 O 2 ' 

at two locations (x , t ) and (x - Ax, t ). 

v o’ o' v o o' 

In a later paper by Rubin and Berstein [19] it was shown for cne- 

phase flow that a modified version of Richtmeyer 1 s 2-step method was 

more stable in solving the system of equations. The use of half-time 

steps was eliminated by averaging F differences at present time, t Q , and 

future time, t + At. This is written as 
o 

Step 1 


,.t+At _ x+Ax * x At ,t „t] 

U x+Ax "2 + Ax x»Ax * r xJ 

2 


(.b.S) 


Step 2 


U 


t+At 

x 


At 

+ 2Ax 


■ 1 

F 

t 1 F C*At 

( ,t*At. ‘ 

2 

x* Ax 

x-Ax 1 x+Ax 

‘ x-Ax 



-1 T 

2 - 


-•9) 


This method is illustrated in Figures D 4 and D.5 on the following page. 

The following is a summary of the schemes described in this Appendix. 
Each will be applied to the gas continuity equation of l-D, two-phase flow. 
Here, T represents the gas being generated from the solid phase. 

Gas continuity equatio n 


at 


ap i u f, 

ax 


(u. 10) 
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